%00000000000000000000000000000000000000000000000000000000000
%THIS FILE SIMULATES A WEEKLY BUFFER STOCK MODEL WITH BI-WEEKLY PAYCHECK
%GENERATES FIGURES 3B,4,6
%00000000000000000000000000000000000000000000000000000000000

%set up environment
    clear all
%print flag (0 for no print, 1 for print)
    print_flag=1;
%set parameters
    R = 1+(0.01/52); %interest rate
    rho = 1; %risk aversion
%three point weekly beta distribution
        beta=[0.91 0.95 0.99];
        n_beta=size(beta,2);
%pay_share (percent of total income)
    pay_share = 0.7;
%temp shock distribution (spread the rest of the payshare over each week)
    %you get 0 on non-paycheck week
    %mu=average uncertain income
    %sigmat - calibrated from previous paper converting monthly to weekly
    %vol
    n_at=5;mu_at=(1-pay_share);sigmat=0.10;
    [theta,p_theta]=tauchen(n_at,mu_at,0,sigmat,4);
    p_theta=p_theta(1,:)';
%paycheck 
    n_p = 2;tp=(1:n_p); % time periods
    mu_at_pay=n_p*pay_share; %paycheck income as a fraction of total
%check that 2 week income adds up to 2
    assert(mu_at_pay + 2*mu_at==2);
%set state space and consumption function
    max_a=15;
    n=1000; %size of grid
    a=(exp(exp(exp(linspace(0.00,log(log(log(max_a+1)+1)+1),n))-1)-1)-1)';  % set up triple exponential grid
%mp(grid points,diff temp shock values * diff paycheck values,week)
    mp = zeros(n,n_at,n_p); %all possible values of next period cash on hand (state x income x pay periods)
%declare functions    
    IUP=inline('c.^(-1/rho)','c','rho'); %incverse MU
    UP=inline('c.^(-rho)','c','rho'); %MU    
%this is the grid of possible m_{t+1} values 
    for i = 1:n_p
        %in the 2nd week get a paycheck so in the first week next week is
        %pay week
        if i == 1
            mp(:,:,i) = a + repmat(theta',n,1)+mu_at_pay;
        else
            mp(:,:,i) = a + repmat(theta',n,1);
        end
    end    
%%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
% Find the steady state consumption function  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
     %define current and next period c
        clear c m;
        c = repmat(a,1,n_p,n_beta); %we start next period c with consuming everything
        c_new = repmat(a,1,n_p,n_beta);
    %define current and next period m because we need it to extrapolate
        m=repmat(a,1,n_p,n_beta);
        m_new=repmat(a,1,n_p,n_beta);      
    %get rid of warning prompts 
    try
        warning('off','last');
    catch
        warning('we cant turn the last warning off yet since it hasnt occured');
    end
 
%interp type
    interp_type='linear';
%loop over different beta vlues     
for b = 1:n_beta         
    diff=1000;
     while diff > 1e-6
        %look over months to pay
        for j = 1:n_p
            %generate index (k) to use for next period consumption function
            if j==n_p
                k=1;
            else
                k = j+1;
            end
        %interpolate the consumption function next period conditional on
        %beta
            if m(1,k,b)==0
                pp = interp1(m(:,k,b),c(:,k,b),interp_type,'pp');
            else
                pp = interp1([0;m(:,k,b)],[0;c(:,k,b)],interp_type,'pp');
            end
        %calculate possible c_{t+1} values based on next period coh    
            cprime = ppval(pp,mp(:,:,j));
        %marginal utility of possible c_{t+1}
            uprime_prime = UP(cprime,rho);
            uprime_prime(uprime_prime>realmax)=realmax;
        %expected marginal utility of c_{t+1}
            e_uprime_prime = uprime_prime * p_theta;
        %back out what c is based on euler equation
            c_new(:,j,b) = IUP(e_uprime_prime * beta(b)* R,rho);
        %calculate cash on hand
          m_new(:,j,b) = a + c_new(:,j,b);
        end
        diff = max(max(abs(c_new(:,:,b)- c(:,:,b)))); %take the largest difference between the two value functions
              c=c_new;m=m_new; 
     end
end
     
   %make sure consumption doesn't exceed coh
        assert(max(max(max(c > m)))==0)   
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GENERATE SIMULATED CONSUMPTION PATH (FIGURE 3B)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%n=people, t=time periods
    discard=100;n = 150;t=1000;T=t*n;discard=discard-1; 
    clear c_sim a_sim y_sim z_sim t_sim tax_refund 
    c_sim=zeros(t,n);a_sim=zeros(t,n);t_sim=zeros(t,n);y_sim=zeros(t,n);z_sim=zeros(t,n);

    %generate random income and betas
        cutoff=cumsum(p_theta',2); %calculate the discrete cutoff points in the distribution of income
        cutoff_beta=[1/n_beta:1/n_beta:1]; %calculate the discrete cutoff points in the distribution of income
        random_y=zeros(T,1);
        rand('seed',1234);
        epsilon=rand(T,1);       
        epsilon_beta=rand(n,1);

     %assigns the random components based on the discretized distribution
      %the idea is to take the N(0,1) draw and see where it falls in the cumulative
      %dist. Then assign a discrete value to it
        idx1=(n_at+1)-sum(repmat(cutoff,T,1) >= repmat(epsilon,1,n_at),2); %look how far we are in the CDF
        idx_beta=(n_beta+1)-sum(repmat(cutoff_beta,n,1) >= repmat(epsilon_beta,1,n_beta),2); %look how far we are in the CDF
            idx_beta(idx_beta>n_beta)=n_beta;
        random_y=reshape(theta(idx1),t,n); %reshape the random shocks back to matrix form
    %create simulated y
        y_sim =  random_y + repmat(ismember(mod(1:t,n_p),0)',1,n).*mu_at_pay; %create the simulated income values
    %start with asset value of 1
            z_sim(1,:)=1;t_sim(1,:)=1;a_sim(1,:)=0;

%generate history of CoH and consumption            
%we have a different consumption function for each period
for b = 1:n_beta
    beta_index = (idx_beta == b);
        for i=2:t
               period_index=mean(t_sim(i-1,beta_index));
            
                %Step 1: What was consumption last period? This is based on the
                %starting value of CoH in period 1
                    pp = interp1([0;m(:,period_index,b)],[0;c(:,period_index,b)],interp_type,'pp');  
                    c_sim(i-1,beta_index)=ppval(pp,z_sim(i-1,beta_index));

                %Step 2: Calculate CoH
                    a_sim(i-1,beta_index) = z_sim(i-1,beta_index) - c_sim(i-1,beta_index);
                    z_sim(i,beta_index) =  y_sim(i,beta_index)  + a_sim(i-1,beta_index);    

                %Step 3: Iterate time period 
                    if t_sim(i-1,beta_index)==n_p
                        t_sim(i,beta_index)=1;
                    else
                        t_sim(i,beta_index)=t_sim(i-1,beta_index)+ 1;
                    end
        end
end

    %discard last observation and the first until a specified point to
    %stabilize the series
        c_sim(end,:)=[];z_sim(end,:)=[];y_sim(end,:)=[];a_sim(end,:)=[];t_sim(end,:)=[];t=t-1;
        c_sim(1:discard,:)=[];z_sim(1:discard,:)=[];y_sim(1:discard,:)=[];a_sim(1:discard,:)=[];t_sim(1:discard,:)=[];t=t-discard;
    %create index for example
        index_short=150:180;
    %calculate average for each beta
         c_time_avg_1=mean(c_sim(:,idx_beta == 1),2);
         c_time_avg_2=mean(c_sim(:,idx_beta == 2),2);  
         c_time_avg_3=mean(c_sim(:,idx_beta == 3),2);  
  %plot time series figure as an example of what the process looks like
  %(figure 3b)
      clf
      hold on
          plot(index_short,log(c_time_avg_1(index_short))-log(mean(c_time_avg_1)),'color','k','LineWidth',2)
          plot(index_short,log(c_time_avg_2(index_short))-log(mean(c_time_avg_2)),'color','k','LineWidth',2,'LineStyle','--')
          plot(index_short,log(c_time_avg_3(index_short))-log(mean(c_time_avg_3)),'color','k','LineWidth',2,'LineStyle',':')
      hold off
      h_legend=legend('\delta=0.91','\delta=0.95','\delta=0.99','Location','southoutside','Orientation','horizontal');
      h_legend.FontSize=12;
      ylabel('Log deviations');
      axis([150,180,-0.1501,0.1501]);
      if print_flag==1 
        print -depsc 'fig\time_series_sim_liq.eps'        
      end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%
%EULER EQUATION TEST [PAY PERIOD] (FIGURE 4)
%%%%%%%%%%%%%%%%%%%%%%%%%%
clear log_avg_liq rel_liq cg cg_pay xx2 ind yy2 cg cg_pay reshape_cg reshape_term reshape_liq
    %calcualte average liquidity conditional on pay or non-pay period
        pay_period_index=(t_sim==2); %index for pay period
        no_payperiod_index=(t_sim==1);        
        log_avg_liq_nopay = log(mean(a_sim(no_payperiod_index(:,1),:))); %avg liquidity for each person in pervious period
            prev_coh = [zeros(1,n);a_sim(1:end-1,:)];
                %make 0 very small numbers
                    prev_coh(prev_coh<=0)=eps;
        rel_liq = log(prev_coh(pay_period_index(:,1),:)) - repmat(log_avg_liq_nopay,size(a_sim(period_index(:,1),:),1),1); %rel liq for each person
                                  
    %calcualte consumption growth
        cg = log(c_sim(2:end,:))-log(c_sim(1:end-1,:)); %consumption grown from t to t+1
        cg_pay=cg(pay_period_index(1:end,1),:); %capture all the pay_periods 

    %set colors and line type
         C = {[0 0 0.5],[0 0.5 0],[0.5 0.5 0.5]}; % Cell array of colors.
         type = {'-','--',':'}; % Cell array of line type.
        
%PLOT THE FIGURE
    clf
    hold on
    for b = [1 2 3]
            beta_index = (idx_beta == b);
            reshape_term = size(rel_liq(:,beta_index),1) * size(rel_liq(:,beta_index),2);
            reshape_liq = reshape(rel_liq(:,beta_index),reshape_term,1);
            reshape_cg = reshape(cg_pay(:,beta_index),reshape_term,1);  

        yy2 = smooth(reshape_liq,reshape_cg,0.90,'lowess');
        [xx2,ind] = sort(reshape_liq);
        plot(xx2,yy2(ind),'color','k','LineStyle',type{b},'LineWidth',2)
    end    
        hold off
        axis([-2,2,-0.21,0.11]) ;
        ylabel('Expected consumption growth');
        xlabel('Relative liquidity');
        h_legend=legend('\delta=0.91','\delta=0.95','\delta=0.99','Location','southoutside','Orientation','horizontal');
        if print_flag==1
            print -depsc 'fig\gc_plot_pay.eps'
        end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%
%EULER EQUATION TEST [NON-PAY PERIOD] (FIGURE 6)
%%%%%%%%%%%%%%%%%%%%%%%%%%
clear log_avg_liq rel_liq cg cg_pay xx2 ind yy2 cg cg_pay reshape_cg reshape_term reshape_liq
    %calcualte average liquidity conditional on pay or non-pay period
        log_avg_liq_pay = log(mean(a_sim(pay_period_index(:,1),:))); %avg liquidity for each person in pervious period
            prev_coh = [zeros(1,n);a_sim(1:end-1,:)];
                 prev_coh(prev_coh<=0)=eps;
            rel_liq = log(prev_coh(no_payperiod_index(:,1),:)) - repmat(log_avg_liq_pay,size(a_sim(period_index(:,1),:),1),1); %rel liq for each person

    %calcualte consumption growth
        cg = log(c_sim(2:end,:))-log(c_sim(1:end-1,:)); %consumption grown from t to t+1
        %add in a row of zeros so the dimensions are correct
            cg=[cg;zeros(1,n)];
        cg_pay=cg(no_payperiod_index(1:end,1),:); %capture all the pay_periods 
 
%PLOT THE FIGURE
    clf
    hold on
    for b = [1 2 3]
            beta_index = (idx_beta == b);
            reshape_term = size(rel_liq(:,beta_index),1) * size(rel_liq(:,beta_index),2);
            reshape_liq = reshape(rel_liq(:,beta_index),reshape_term,1);
            reshape_cg = reshape(cg_pay(:,beta_index),reshape_term,1);


        yy2 = smooth(reshape_liq,reshape_cg,0.90,'lowess');
        [xx2,ind] = sort(reshape_liq);
        plot(xx2,yy2(ind),'color','k','LineStyle',type{b},'LineWidth',2)

    end    
        ylabel('Expected consumption growth');
        xlabel('Relative liquidity');
        h_legend=legend('\delta=0.91','\delta=0.95','\delta=0.99','Location','southoutside','Orientation','horizontal');
        if print_flag==1
            print -depsc 'fig\gc_plot_nopay.eps'        
        end